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Abstract. The J'-'-statistic is an optimal detection statistic for continuous 
gravitational waves, i.e., long-duration (quasi-)monochromatic signals with 
slowly-varying intrinsic frequency. This method was originally developed in the 
context of ground-based detectors, but it is equally applicable to LISA where many 
signals fall into this class of signals. We report on the application of a LIGO / GEO 
.7-"-statistic code to LISA data-analysis using the long- wavelength limit (LWL), and 
we present results of our search for white-dwarf binary signals in the first Mock 
LISA Data Challenge. Somewhat surprisingly, the LWL is found to be sufficient - 
even at high frequencies - for detection of signals and their accurate localization 
on the sky and in frequency, while a more accurate modelling of the TDI response 
only seems necessary to correctly estimate the four amplitude parameters. 
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1. Introduction 

The Mock LISA Data Challenge (MLDC) 1 has the purpose of encouraging the 
development of LISA data-analysis tools and assessing the technical readiness of 
the community to perform gravitational-wave (GW) astronomy with LISA. The first 
round of the MLDC was released in June 2006 |2j, the submission deadline was in 
December 2006 and a report summarizing the submitted results has been published [3] . 
The challenges consisted of several data-sets containing different types of simulated 
sources and LISA noise. The three types of sources are white-dwarf binary signals 
(WD), coalescing supermassive black holes (SMBHs) and extreme mass-ratio inspirals 
(EMRIs). 

The data analysis of LISA poses a few specific difficulties not encountered 
in ground-based detectors: the signal (reduced) wavelength is typically not long 
compared to the arm-length of the detector, so the long-wavelength limit (LWL) 
does not generally apply. Furthermore, in order to cancel the dominating laser- 
frequency noise, one has to analyze intricate algebraic combinations of time-delays 
between spacecraft instead of simple "strain" , an approach known as time-delay 
interferometry (TDI). Another difficulty stems from the large number of detectable 
sources in the LISA bandwidth, which complicates their separate detection and 
parameter estimation, usually referred to as the "confusion problem" . 

Most of the relevant signals in LISA (WD, SMBH, EMRI) will be long-lasting 
(of the order of a year) and are (quasi-)monochromatic with slowly-varying intrinsic 
frequency /(r); in this sense they belong to the class of continuous GWs. In the case 
of ground-based detectors the typical sources of continuous GWs are spinning neutron 
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stars with non-axisymmetric deformations. One of the standard tools developed for 
these searches is the J-"-statistic |3], which is an optimal detection statistic (in the 
sense of the Neyman-Pearson lemma) based on matched filtering. We have restricted 
our searches in the first MLDC to WD-binary signals, which are very similar to GWs 
from spinning neutron stars, which have very little intrinsic frequency evolution / 
(in fact, here it was / = 0) and constant orientation and polarization. Contrary to 
the approach used in O [6] , we use an .F-statistic code developed for the continuous- 
wave search in LIGO/GEO, with only minimal modifications to adapt it to LISA. In 
particular, we use the LWL at all frequencies, which turns out to work surprisingly 
well even at high frequencies where the reduced wavelength is comparable to the LISA 
arm length. 

2. Methods and Pipeline 

2.1. Continuous Gravitational Wave Signals 

A system with an oscillating mass quadrupole moment emits GWs described, far from 
the source, by the metric perturbation 

h^ A+ cos{(j}Q + (f>)*e++ Ay, sin(0o + 0) e^x , (1) 

where ^+ = ex®ex — ey®ey and e*x = Cx^ey+ey^Cx are the polarization basis tensors 
constructed from a right-handed basis {ex, ey, ez} with e^ pointing in the direction 
of propagation of the wave, described by the ecliptic latitude P and longitude A, and 
Cx and By along the principal polarization axes. In an inertial reference frame, such 
as the solar-system barycenter (SSB), the phase of this (quasi-)periodic signal can 
be written as 0(t) — 2tt ^ fi^') dr' , in terms of the (slowly-varying) intrinsic GW 

frequency /(r) = /(rj-ef ) + /(TVef ) Ar -|- . . . , where Trcf is a reference time at which the 
frequency and spindown parameters are defined, and At = r — Trcf. The WD signals 
in the first MLDC were restricted to have a constant intrinsic frequency, i.e., /(r) = /. 
This is a realistic assumption at low frequencies / ~ 1 mHz, but at higher frequencies 
/ ^ 10 mHz one would probably have to include one derivative / (e.g. see j5j) in an 
actual search on LISA data. In the case of a binary system for which orbital evolution 
due to GW emission can be neglected, the principal polarization axes are found by 
taking the unit vector ex to lie in the orbital plane and ey in the hemisphere containing 
the orbital angular momentum. The polarization amplitudes are = /io(l-(-cos^ t)/2 
and Ax = ho cos where ho is usually referred to as the amplitude of the GW, and 
i is the inclination angle between the propagation direction ez and the normal to the 
orbital plane. In order to separate the sky position {(3, A} from the source polarization, 
it is useful to consider a polarization basis associated only with the sky position; 
this is done by defining a right-handed orthonormal basis {e^, e,,, 6,^} with e,^ — ez 
as the propagation direction, lying in the ecliptic plane and e^ in the northern 
hemisphere. The alternative polarization basis is then £*+ = (g) — (8" and 
£x = e(^® ejf-V e^ (g) Cit, and the principal polarization axes of the GW are determined 
by the angle -0 from to e^, measured counter-clockwise around e-z = ec, i-e., 

e*+ — £*+ cos 20 -I- e^x sin 20, ^ 
e*x = — £*+ sin 20 -|- e*x cos 20 . 
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Figure 1. LISA configuration and TDI conventions used. 

In terms of this alternative polarization basis, the GW tensor can be written as 

4 

Mr)^^^''/:^(r), (3) 

where the four amplitude parameters {A'^} are 

A^ = cos (po cos 2ijj — Ax sin (po sin 2ip , 

A^ — cos 00 sin 2^ + Ax sin cos 2iIj , 

^ (4) 
A = —A+ sin (po cos 2ijj — Ax cos (/)o sin 2-0 , 

= — j4+ sin 00 sin 20 + Ax cos 0o cos 2-0 , 

while the tensors {/i^} depend on the frequency /(t) and the sky position A}, 
namely 

hi{T) = 7+ cos0(t) , /i2(t) = e'x cos0(r) , 

^ ^ ^ ^ (5) 

/i3(r) = e+ sin0(r) , ft,4(T) = ex sin0(r) . 

Note that the geometrical conventions for the amplitude parameters used here are 
consistent with the LIGO/GEO conventions for continuous GWs (e.g. 7 ), but differ 
from the LISA/MLDC conventions [31 [8]. The translation into MLDC conventions is 
given by: "Amplitude" = ho/2, u ^ -k — u, ip ^ 7r/2 — ip and 0o ^ 0o + tt. 

2.2. LISA Response in the Long-Wavelength Limit 

The LISA design consists of three spacecraft with laser links between each pair, in 
a geometry illustrated in figure [TJ The MLDC data were generated by two different 
programs: Synthetic LISA [9] simulates a detector output consisting of Doppler shifts 
of the LISA lasers due to relative motion of the spacecraft, while LISA Simulator [10] 
simulates the phase differences between laser light following different paths between the 
spacecraft. In both cases the underlying variables are combined with appropriate time 
shifts to form TDI observables which cancel the (otherwise dominating) laser frequency 
noise [iTJ O [5] . One choice of such TDI quantities is the set of three observables 
{X, Y, Z}, which were used to pubHsh the data of the first MLDC. These observables, 
which can be thought of as representing the output of three virtual "detectors" 
J, are related to the gravitational wave h through somewhat involved expressions 
depending on the frequency and propagation direction of the wave. However, in the 
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LWL approximation, in which the reduced wavelength c/(27r/) is assumed to be large 
compared to the distance between the spacecraft, i.e., / <C lOmHz, the responses can 
be approximated (assuming Li k L2 ~ k, L) as 

^^synthLISA = _ . ^ ^LISAsim :h, (6) 

c 

where : denotes the contraction of both tensor indices, and = {n2 ®n2~n^® ns) /2 
is the usual LWL response tensor for a GW interferometer with arms n2 and 713. The 
analogous expressions for Y and Z are obtained by cyclic permutations of the indices 
1 — > 2 — > 3 — > 1. We define an associated scalar "strain" for each of the detectors 

I ^ X,Y,Z as 

h\t)^l\t):h{T{t)) . (7) 

The timing relation r(t) accounts for the Doppler effect caused by the orbital motion 
of the detector, namely r(t) = t — f- Sz/c, where r{t) is the position of the detector 
with respect to the SSB, and Cz is the propagation direction of the GW. Note that 
in the LWL approximation, we can assume that all virtual detectors follow the same 
trajectory r{t) corresponding to the barycenter of the three spacecraft. 

The input to our search code consists of Fourier-transformed data stretches 
of duration Tsft, the so-called SFTs (for Short Fourier Transforms), which is a 
common data format used within the LIGO Scientific Collaboration for continuous- 
wave searches (e.g. see [7]). The time baseline Tsft has to be chosen sufficiently short 
such that the noise-floor can be approximated as stationary and the rotation and 
acceleration of the detector can be neglected. For LISA we chose Tsft = 7 days, while 
in LIGO / GEO (where the rotation of the Earth dominates the acceleration) this is 
typically Tsft = 30min. Approximating the detector tensor as constant during 
TsFT) we can Fourier-transform ([6]) to obtain 

We use h^ (/) as our SET input data, and so we can run the same pipeline on data from 
LISA Simulator and synthetic LISA, with only a different "calibration" ([5]) used to 
generate the SFTs. The noise contributions to X, Y, and Z are correlated, therefore 
it is often convenient to work with the TDI variables X and Y — Z instead, which are 
statistically independent. This is a straightforward generalization, using the response 
tensor d^~^ = — d^ . Note that X and Y— Z generally have different noise levels, 
but this is properly taken into account in the multi-detector J^-statistic. 



2.3. The T -Statistic Method 

The JF-statistic was originally developed in [4] and extended to the multi-detector case 
in [T3j. A generalization to the full TDI framework for LISA was obtained in but 
here we follow the simpler route of working in the LWL approximation, which allows 
for a more direct application of existing LIGO/GEO codes to LISA data analysis. 

Combining the scalar strain ([7]) with the expression ([3]) for the GW tensor, we 
can write the strain signal at detector / as 

h\t)^j2-^'K(*^' (9) 
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in terms of the four basis functions 

=a^(i) cos0(r(i)), hi{t)^h\t) cos ct){T{t)), 
hl{t) ^a\t) sin (I) {T{t)), hiit) = b\t) sin (j){T{t)) , 

where we defined the antenna-pattern functions = : e+ and = : e'x- 
The functions {/i^} depend on the sky-position {f3, A} and the frequency /(r) of 
the source. We see that the signal parameters separate into two classes: (i) the 
four amplitude parameters A = {A'^} given in ([4]) and (ii) the Doppler parameters 
9 = A, /,/,/, ...}. We model the output {t) of detector / as a superposition 
of stationary Gaussian noise (t) and a signal {t; A, 9) . Following the notation of 
[131 [Uj, we write the different data-streams {t) as a vector a;(t), and we define the 
standard multi-detector scalar product as 

/oo 
^'*{I)Si}{f)y'{f)df, (fl) 
i,j 

where x is the Fourier-transform, x* denotes complex conjugation, and {<S'7j (/)} are 
the elements of the inverse of the noise-power matrix. We search for a signal by seeking 
the parameters {A, 0} which maximize the log-likelihood ratio 

L{x- A 9) = {x\h) - ]^{h\h) = A^'{x\h^:) - ]^A^'{h^,\K)A'' , (12) 

with automatic summation over repeated amplitude indices /i, v. Defining 

x^,{9) = {x\h^), and = , (13) 

we see that L is maximized for given 9 by the amplitude estimator A'^^^^^ — M'^'^'^Xi,, 
where A4'^'^ is the inverse matrix of A^p,y. Thus the detection statistic L, maximized 
over the amplitude parameters A, is 

T{x;9) = ^x^M^''x,, (14) 

which defines the (multi-detector) J^-statistic. 



2.4- Analysis Pipeline 

Our analysis was based on standard LAL/LALApps software [T5j developed for 
the search for continuous GWs with ground-based detectors, in particular the code 
ComputeFStatistic_v2, which implements the multi-detector J^-statistic ()14p . Only 
minor modifications were necessary to adapt this code to the analysis of LISA data 
using the LWL approximation. All white-dwarf binary signals in the first MLDC 
had constant intrinsic frequency /, so the set of Doppler parameters to search over 
consisted of 9 — {/3, A,/}. We performed a hierarchical search that first runs single- 
detector searches on each of the TDI variables /, looks for coincident local maxima of 
2J^, and in the last step performs a multi-detector JF-statistic search to establish the 
parameters of each candidate signal. Our initial analysis submitted as an MLDC 
entry [3] used the TDTvariables X, Y and Z as three "detectors", assuming for 
simplicity that their correlation matrix Sij{f) is diagonal. However, given that 
the corresponding noises are correlated, we subsequently re-ran the search using the 
uncorrelated TDI variables X and 1" — Z, which was used for the results presented 
here (but did not result in any significant changes in the results). Whether / ranges 
through {X, Y, Z} or {X, Y — Z}, the structure of the pipeline is the same: 
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(i) Perform a wide-parameter J^-statistic search on each data stream / over a 
template grid of Doppler parameters A, /} . The grid was chosen as isotropic 
in the sky, with angular mesh size da = \/2m/(27r/i?oib/c), with the orbital 
radius i?orb = 1 AU and we a mismatch of m = 0.3. The frequency spacing used 
is df = Vl2 ml {t: T), where T — 1 y is the observation time. These step sizes 
were computed from the orbital metric |T4] . 

(ii) Keep only candidates which are local maxima of 2J- (above some threshold), and 
which are coincident with consistent Doppler parameters in all detectors /. 

(iii) Perform a more finely-gridded multi-detector search around each candidate to 
increase the accuracy of the parameter estimation. 

(iv) Classify each candidate as primary if it has the highest 2J- value within 
A/ ~ 1.4 X 10^^/, and as secondary otherwise. 

The last step arises from the empirical observation that a given signal will have 
secondary "false" JF-statistic maxima at frequencies within roughly ~ 10~^ / but at 
different sky positions. Only primary candidates were reported, while the secondary 
candidates were discarded. This is a limitation of our pipeline: given two signals 
very close in frequency but at different sky positions, it cannot distinguish the peak 
at the true sky position of the "fainter" source from a secondary maximum of the 
"brighter" one. This problem is seen particularly in Challenges 1.1.4 and 1.1.5 with 
signals clustered very densely in frequency. 

3. Results 

3.1. Challenge 1.1.1: Isolated Binaries 

This challenge consisted of three separate data sets, each containing one WD signal at 
an unspecified sky position and within a given frequency band: in 1.1.1a at ^ 1 mHz, 
in 1.1.1b at ^ 3 mHz, and in 1.1.1c at ~ 10 mHz. Note that the LWL is only 
a good approximation for / <C 10 mHz, and we therefore expect it to deteriorate 
significantly in 1.1.1b and 1.1.1c. Nevertheless, in each of the three cases our pipeline 

Table 1. Recovery of Doppler parameters in Challenge 1.1.1: A/ is the 

frequency error, and (/igky is the angle between recovered and true sky position. 
Frequency / and sky position {/3, A} were accurately determined even at the 
highest frequencies. 



Challenge 


/ 


13 


A 


A/ 


"/■sky 


1.1.1a 
1.1.1b 
1.1.1c 


1.1 mHz 
3.0 mHz 
10.6 mHz 


0.95 rad 
-0.09 rad 
-0.11 rad 


5.07 rad 
4.63 rad 
4.66 rad 


1.7nHz 
0.8 nHz 
0.2 nHz 


34.8 mrad 
7.1 mrad 
4.4 mrad 



recovered a single primary candidate, and the Doppler parameters were determined 
with very good accuracy, as summarized in table [H The apparent improvement in the 
Doppler accuracy seen in table [1] is due to statistical fluctuations. Running this search 
on a larger number of sources (such as in Challenge 1.1.2) reveals no clear trends. 
The recovery of the amplitude parameters A is illustrated in figure [2l comparing 
the estimated 4-vector ^cand to the 4-vector ^kcy of the injected parameters. The 
amplitude 4- vectors A live in a space with constant metric tensor Ai^y given in (fT^l) . 
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Figure 2. Recovery of amplitude parameters in Challenges 1.1.1a {left), 1.1.1b 
{middle), and 1.1.1c {right). Each plot compares the recovered amplitude 4- vector 
-4cand to the injected signal 4- vector ylkoy i shown in the plane defined by the two 
vectors. Gaussian fluctuations would lead to a separation of the endpoints of the 
order | Ayl| ~ 2. The breakdown of the LWL with increasing frequency leads to 
larger errors, affecting both the orientation and the magnitude of the recovered 
amplitude vector. 



SO the norm is \A\^ = A'-' M^u A^. The two vectors .4cand and .4key define a plane, and 
so we can plot them in two dimensions, with the horizontal and vertical components 

•^kcy ' A 
l^kcy I 



^11 = , / , , Ai 



A - A, '^'^"^ 



(15) 



l-^kcy I 

where the inner product is calculated using the metric AAp^u- These components 
are shown in figure [U and we see that the agreement of the amplitude parameters 
deteriorates substantially for higher frequencies, where the LWL approximation breaks 
down. 

If the deviation is caused by noise alone, then the difference = ^cand ~ ^kcy 
between the amplitude vectors has zero mean, i.e., E [i^A\ = 0, and covariance 
£^[A^^ A^'^] = M.^^ , where E[. . .] denotes the expectation value. The magnitude 
of this difference, \l\A\ = ^Ay4,^ A^^^ I^A" , would have variance £'[|A^|^] = 
M.^iuJ^^'^ = 4. Therefore |A^| /2 measures the difference between the two amplitude 
vectors in terms of the number of standard deviations. It is also instructive to compare 
the magnitude of the recovered versus the injected amplitude vector. The magnitude 
l^kcyl of the injected signal is equivalent to the optimal signal-to-noise ratio (SNR). 
Note, however, that |^cand|^ is a biased estimator for |^kcy|^, namely 



E 



lAandl' = E [2J^] = 4 + I Acyr ■ (16) 



|2 



Therefore we use the following measure for the error in the norm of the recovered 
amplitude vector -4cand: 

_ |-4cand| ^ l-^kcyl ~4 

— — ■ ™ 

which is unbiased, i.e., [e^] = 0. The standard deviation of |ylcand| is 2(2 + 
|-4key| ~ 2|Acy|, and so the expected error from noise alone would be 
■£'[£.4] ~ l-^keyl ^- Table [5] summarizes the errors in the amplitude parameters for 
the three challenge data sets in terms of \AA\ /2, the relative difference of the 
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Table 2. Errors in recovered amplitude parameters in Challenge 1.1.1: as seen in 
figure [2] the angle i/i_4 between ylcand a^nd -^key grows with increasing frequency, 
and there is an increasing deficit in the magnitude l-A^andl with respect to the 
SNR |^kcy|i ^ quantified by tj^. The absolute error in |A^| /2 from Gaussian 
noise would be expected to be ~ 0(1), while for and (f)j^ it would be |.4i;cy[ ^ ■ 



Challenge 


/ 


Mkcyl"' 




4>A 


|A^| /2 


1.1.1a 
1.1.1b 
1.1.1c 


1.1 mHz 
3.0 mHz 
10.6 mHz 


0.020 
0.027 
0.007 


0.005 
-0.020 
-0.306 


0.059 
0.527 
2.207 


1.5 
9.6 

108.7 



norms, and the angle <j)j( between the recovered and the injected amphtudc vectors, 
given by 

We see in table [5] that the amplitude errors are larger than would be expected from 
noise fluctuations alone, especially at higher frequencies, which is consistent with the 
breakdown of the LWL. 



3.2. Challenge 1.1.2: Verification Binaries 

In Challenge 1.1.2, the sky position and frequency of twenty "verification binaries" 
was given, while the amplitude parameters of the injected signals were unknown. 
We therefore performed a targeted J^-statistic search at each of the specified sets 
of Doppler parameters, and found the maximum-likelihood estimators ^cand for the 
amplitude parameters. Figure [3] illustrates the discrepancies between the recovered 
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Challenge 1.1.2 Parameter Errors 
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Figure 3. Recovery of amplitude parameters in Challenge 1.1.2. Left: all 20 
signals are recovered with l^candl ^ 8.5, but the errors |A^| /2 are substantially 
larger than the expected standard deviation of unity for all but the smallest 
frequencies. Top right: the angle 4>A between the true and recovered amplitude 
vectors grows with frequency, and is always larger than its expected standard 
deviation of jytkoyl ■ Bottom right: the norm of the recovered amplitude vector 

is within the expected range of |.4]jey| ^ for much of the frequency band, but 
begins to show a deficit for / > 5 mHz. 



-4cand and the injected amplitude parameters Aoy, in terms of |Ayl|, e^, and 4>a- 
Again we see that our recovered amplitude parameters differ from the injected ones 
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by more than would be expected from Gaussian noise alone, and that the agreement 
deteriorates at higher frequencies. 



3.3. Challenge 1.1.3: Resolvable Binaries 

Challenge 1.1.3 was a blind search on data containing 20 white dwarf binary signals 
across the LISA band. As shown in figure HI we recovered 17 of the 20 signals with 
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Figure 4. Dopplcr parameter recovery in Challenge 1.1.3: errors in sky position 
{left figure) and frequency {right figure) as functions of frequency. The width of 
the frequency coincidence window, i.e., 1.4 X 10~*/, is too small to be seen on 
this scale, so no error bars are shown on the found injections. The three missed 
signals (long vertical lines) all fall close to recovered signals, but outside of all 
coincidence windows. 



good frequency and sky accuracy. The three missed signals were at frequencies "close" 
to recovered sources, but not within the frequency coincidence window of 1.4 x IQ^^f, 
and there is some indication that the Doppler parameters of those sources were slightly 
compromised. 



3.4^. Challenge 1.1. 4 and 1.1.5: Source Confusion 

In Challenges 1.1.4 and 1.1.5, many sources were injected into a small frequency 
band in order to illustrate the source confusion problem, namely 45 signals within 
[3, 3.015] mHz in Challenge 1.1.4 and 33 signals within [2.9985, 3.0015] mHz in 
Challenge 1.1.5. As shown in Figure El our pipeline "found" signals all across the 
band, namely 25 signals in Challenge 1.1.4 and only 5 signals in Challenge 1.1.5, but 
many of them were far removed in sky position from any true signal. Many additional 
signals were missed within the frequency coincidence window, presumably because 
they were mistaken for secondary maxima of the "found" signals. The results of this 
challenge illustrate a known limitation of the pipeline used here: it cannot distinguish 
multiple signals too close together in frequency. 



4. Conclusions 



Using the J^-statistic in the long-wavelength limit approximation, we found that 
the estimation of the four amplitude parameters {A'^} deteriorates significantly with 
increasing frequency, as would be expected from the breakdown of the LWL. However, 
the detection of signals and the estimation of the Doppler parameters (frequency and 
sky-position) does not seem to be affected by the use of the LWL, even at frequencies 
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Figure 5. Doppler parameter recovery in Challenges 1.1.4 {Top row) and 1.1.5 
{bottom row): errors in sky position {left column) and frequency {right column) as 
functions of frequency. Each of the "missed" injections falls within the coincidence 
window of a recovered signal, and would therefore have been rejected as a 
secondary maximum. The dashed lines in the top-right plot show the maximum 
possible frequency recovery error, namely the width of the coincidence window. In 
both challenges, source confusion causes our pipeline to find a candidate at every 
possible frequency, including one false alarm at / = 3.0022mHz in Challenge 1.1.4. 



as high as f lOmHz. This somewhat surprising result suggests the fohowing 
"hierarchical" search strategy: start with a fast JF-statistic code using the LWL to 
detect signals and localize them in Doppler space, then use a more accurate (and 
computationally expensive) modelling of the TDI responses to estimate the amplitude 
parameters. 

We are planning to study these findings in a more systematic way using larger 
number of signals. More work is required to deal with "source confusion", i.e., signals 
that lie within a frequency window 0(1O~^ /). Secondary maxima in parameter space 
due to a signal cannot easily be distinguished from primary maxima corresponding 
to other signals within this frequency window. One popular strategy consists of 
successively "removing" detected signals from the data, which also eliminates its 
associated secondary maxima, and allows one to re-run the search for the next-loudest 
candidates. An alternative approach might consist of a classification of candidates 
into equivalence classes consistent with the same signal, either by using the metric or 
a suitable global correlation criterion analogous to the "circles in the sky" [16] present 
for short observation times. 
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